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1.  Objective 


Boron  carbide  (BC),  due  to  its  extreme  hardness  and  low  density,  has  been  used  as  an  armor 
ceramic  for  many  years.  Its  elastic  properties,  which  surpass  more  dense  compounds  such  as 
silicon  carbide  by  a  factor  of  2,  suggest  that  BC  should  sustain  high  impact  pressures  without 
failure;  however,  there  is  an  anomalous  loss  of  impact  resistance  within  the  material,  apparently 
due  to  increased  fragmentation  by  a  factor  of  2  to  3.  Recent  experimental  work  within  the  U.S. 
Anny  Research  Laboratory  has  identified  the  formation  of  nanoscale-sized  intragranular 
amorphous  bands  as  the  primary  damage  mechanism  leading  to  the  marked  reduction  in  ballistic 
perfonnance  of  BC  (1 ,  2).  The  Army  has  a  need  to  understand  the  mechanism  for  the  fonnation 
of  these  amorphous  bands,  at  the  atomic  level,  which  will  allow  for  the  design  of  chemically 
modified  BC  materials  that  will  show  improved  impact  resistance  and,  hence,  improved  armor 
ceramic  perfonnance.  Using  first  principles  quantum  mechanics,  we  provide  much  needed 
insight  into  the  formation  mechanism  of  the  experimentally  observed  amorphous  bands  that  are 
known  to  weaken  BC  ceramics.  In  particular,  the  stress-induced  amorphization  of  BC  is 
explored  via  application  of  the  -Bom  instability”  criterion,  which  requires  that  the  matrix  of 
second  derivatives  of  the  energy  per  unit  cell  be  positive  definite,  where  the  six  independent 
variables  are  strain  measures  in  the  current  configuration.  A  material  that  exhibits  a  Bom 
instability  loses  its  ability  (at  some  stress  and  temperature)  to  convert  certain  incremental  strain 
patterns  into  stress  increments.  These  strain  patterns  are  known  as  — sdf  modes  of  deformation, 
which  must  be  related  to  movements  within  the  atomic  structure.  Identification  of  the  onset  and 
nature  of  the  Bom  instability  will  give  further  insight  into  the  fonnation  mechanism  of  the 
amorphous  regions  that  weaken  the  BC  material.  In  this  research,  we  examine  the  influence  of 
both  pressure  and  shear  stress  on  the  Bom  instability. 


2.  Approach 


BC,  with  nominal  stoichiometry  B4C,  consists  of  12-atom  icosahedra  crosslinked  by  3-atom 
chains  (figure  1)  and  belongs  to  the  crystallographic  space  group  R-3m  (3).  The  exact  nature  of 
the  3 -atom  chain  is  unknown  experimentally  due  to  the  close  spectroscopic  signatures  of  boron 
and  carbon.  The  crystal  structure  consists  of  a  mixture  of  B12  icosahedra  linked  by  C-C-C  chains 
and  BnC  icosahedra  linked  with  C-B-C  chains  (among  others).  For  the  C-B-C  -polytype,”  there 
is  an  additional  complication;  the  location  of  the  carbon  atom  in  the  BnC  cage  is  not  known. 
This,  again,  is  due  to  the  difficulty  in  definitively  resolving  boron  and  carbon;  therefore,  the 
carbon  atom  inserted  into  the  cage  could  reside  in  one  of  two  symmetry  unique  centers,  the  so- 
called  polar  and  equatorial  sites  within  the  icosahedron,  as  shown  in  figure  1 .  It  should  be  noted 
that  although  BC  is  generally  regarded  to  have  R-3m  symmetry,  this  is  only  tme  for  the  C-C-C 
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polytype  and  placement  of  a  carbon  atom  within  the  icosahedra  causes  a  distortion  of  the 
rhombohedral  lattice  and  the  actual  structure  is  dictated  by  an  averaging  of  all  the  different 
polytypes. 


Figure  1.  Icosahedral  structure  of  BC.  Polar  and 

equatorial  sites  in  the  icosahedra  are  indicated 
by  dark  green  and  white  spheres,  respectively. 

As  a  route  to  the  elucidation  of  the  mechanism  by  which  pressure  and  shear  induce  an 
amorphous  transition  in  BC,  we  have  applied  the  Born-stability  criterion  to  each  of  the  three 
polytypes  discussed  earlier.  Born  showed  that  an  expansion  of  the  internal  energy  of  a  crystal  in 
a  power  series  in  the  strain,  along  with  the  imposition  of  positivity  of  the  energy,  leads  to 
restrictions  on  the  relative  magnitudes  of  the  elastic  constants  of  a  stable  crystal  (4).  The 
existence  of  such  a  Born  instability  can  be  determined  quantum  mechanically. 

BC  is  highly  anisotropic  elastically  and  for  crystals  with  R-3m  symmetry,  there  are  only  six  non¬ 
zero  independent  elastic  constants  { Cq }  for  the  unloaded  material  (5).  These  constants  have 
been  determined  experimentally  for  BC,  except  for  the  modulus  C14  (6).  Each  of  the  elastic 
moduli  varies  independently  with  pressure,  and  at  some  point,  the  system  may  reach  a  structural 
instability.  Imposition  of  the  Bom  stability  conditions  leads  to  the  following  restrictions  on  the 
elastic  constants  for  a  stable  crystal  with  R-3m  symmetry: 


Cn-|C12|>0 

(1) 

(Cn  +  C12)C33  -  2C13*C13  >  0 

(2) 

(Cn  -  C12)C44  -  2C14*C14>  0 

(3) 

The  procedure  previously  adopted  by  others  is  to  compute  the  six  elastic  constants  above  as  a 
function  of  pressure  and  evaluate  equations  1  through  3  to  determine  the  onset  of  the  instability, 
i.e.,  the  pressure  at  which  one  or  more  of  the  equations  above  are  no  longer  valid  (7,  8).  Once 
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the  initial  instability  has  been  located,  evaluation  of  the  -soft  modes”  of  deformation  (atomic 
displacements  corresponding  to  the  instability)  can  be  determined. 

Equations  1-3  are  generally  not  applicable  to  the  incremental  moduli  of  crystals  under  finite 
loading;  however,  a  general  indicator  of  an  elastic  instability  for  all  materials  under  any  state  of 
loading  is 

det  |  C  |  =  0,  (4) 

where  C  is  the  full  6x6  tensor  of  incremental  elastic  moduli  (9,  10).  Further,  since  the 
detenninant  of  a  matrix  is  equal  to  the  product  of  its  eigenvalues,  an  equivalent  statement  of 
equation  4  is  that  all  the  eigenvalues  of  C  must  be  positive  and  the  first  stress  state  at  which  an 
eigenvalue  equals  zero  corresponds  to  an  instability.  It  should  be  noted  that  for  a  crystal  under 
load,  the  elastic  constant  tensor  C  should  be  replaced  by  the  -effective”  elastic  constant  tensor  B 
given  by 

Bijki  =  Qjki  +  (l/2)(6ikiji  +  6jkiii  +  SiiXjk  +  §jiXik  -26kiXij)  (5) 

with  Xij  being  an  element  of  the  stress  tensor  (10). 

The  elastic  constants  were  computed  using  the  CP2K  software  package  (11).  Density  functional 
theory  using  the  Perdew-Burke-Ernzerhof  (PBE)  functional  (12)  in  a  double  zeta  valence  plus 
polarization  basis  set  was  used  for  all  calculations  with  a  planewave  cutoff  of  800  Rydberg. 
Elastic  constants  are  related  to  the  second  derivative  of  the  total  energy  with  respect  to  strain,  s;, 
via 


1  d2E 

V  dsidsi 
'  1  o 


(6) 


where  V  is  the  unit  cell  volume  and  i,j=l . .  .6  in  the  Voigt  notation  (13).  The  CP2K  software 
does  not  evaluate  elastic  constants  analytically;  therefore,  a  Fortran  program  was  written  for  this 
work  that  evaluates  the  second  derivative  of  the  energy  with  respect  to  strain  (equation  6)  via  a 
finite  difference  of  analytic  first  derivatives  provided  by  the  CP2K  code.  This  new  program 
includes  the  required  stress  corrections  from  equation  5  in  the  computation  of  the  elastic 
constants,  computes  the  bulk  modulus  using  the  Voigt-Reuss-Hill  approximation  (14),  and 
finally  detennines  the  eigenspectrum  of  the  elastic  tensor.  Since  we  have  used  the  finite 
difference  approach,  th q  full  6x6  elastic  constant  tensor  can  be  evaluated  with  only  6  quantum 
mechanical  (QM)  calculations  (12  if  double-sided  differences  are  taken,  as  is  done  in  this  work). 
This  is  much  less  than  the  minimum  43  calculations  required  using  a  strain  energy  approach 
advocated  by  some  authors. 
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3.  Results 


The  optimized  lattice  parameters  for  each  polytype  are  presented  in  table  1 .  The  theoretical 
structures  computed  using  PBE  are  in  excellent  agreement  with  experiment  ( 3 )  and  the  distortion 
of  the  lattice  from  purely  rhombohedral  symmetry  is  clearly  evident  in  the  polar  and  equatorial 
polytypes.  Lazzari  et  al.  conducted  a  computational  study  and  concluded  that  the  polar  polytype 
is  the  most  stable  configuration  energetically  (15).  This  is  also  supported  by  our  results  where  it 
is  found  that  the  C-C-C  and  equatorial  polytypes  lie  71  and  35  meV/atom  higher  in  energy  than 
the  polar  configuration,  respectively. 

Table  1.  Unit  cell  parameters  of  BC  polytypes  (lengths  in  angstrom,  angles  in  degrees  volume 
in  cubic  angstroms). 


Polytype 

a 

b 

c 

alpha 

beta 

gamma 

Volume 

ccc 

5.196 

5.196 

5.196 

66.00 

66.00 

66.00 

112.12 

Polar 

5.070 

5.215 

5.215 

65.24 

66.07 

66.07 

109.75 

Equatorial 

5.176 

5.213 

5.213 

64.87 

64.96 

64.96 

110.17 

Experiment 

5.19 

5.19 

5.19 

65.18 

65.18 

65.18 

110.02 

Each  polytype  was  subjected  to  hydrostatic  compression  and  the  resulting  pressure -volume  data 
was  fitted  to  the  third  order  Birch-Mumaghan  equation  of  state  (16).  As  shown  in  figure  2,  the 
pressure  response  of  each  polytype  is  very  similar.  The  computed  bulk  moduli  for  each  polytype 
(and  its  associated  pressure  derivative  resulting  from  the  Birch-Murnaghan  fit)  are  presented  in 
table  2  and  compare  exceptionally  well  with  the  experimental  value.  The  polar  and  equatorial 
polytypes,  both  with  C-B-C  chains,  have  bulk  moduli  that  are  marginally  different;  however,  the 
C-C-C  polytype  has  a  reduction  of  ~12  GPa,  indicating  that  it  is  softer  than  the  other  polytypes. 

It  has  been  suggested  that  the  C-C-C  polytype,  a  minority  phase  in  the  BC  crystal,  is  the 
-weakest  link”  in  the  structure  and  may  be  the  polytype  that  fails  first  when  impacted.  The  lower 
bulk  modulus  is  in  support  of  that  hypothesis  however  more  rigorous  evidence  will  be  provided 
in  the  following. 
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Table  2.  Bulk  modulus  and  pressure  derivative  resulting  from  Birch -Murnaghan  equation  of 
state  fit  to  DFT  data. 


Polytype 

Bulk  Modulus 

(GPa) 

Pressure  Derivative 

ccc 

222 

3.5 

Polar 

234 

3.3 

Equatorial 

233 

3.4 

Experiment 

236 

The  elastic  moduli  of  BC  have  been  evaluated  experimentally  by  McClellan  et  al.  (6);  however, 
in  their  analysis,  they  assumed  a  hexagonal  symmetry  (due  to  limitations  in  their  fitting 
software),  which  leaves  the  value  of  the  C14  elastic  constant  indeterminate  (5).  Since  elastic 
constants  are  tensor  quantities  that  depend  on  orientation,  in  order  to  have  a  direct  comparison 
with  experiment,  we  have  used  a  hexagonal  setting  of  the  rhombohedral  BC  unit  cell  in  the  initial 
QM  calculations.  This  does  not  change  any  of  the  conclusions  resulting  from  the  stability 
analysis  since  the  eigenvalues  are  invariant  to  any  orthogonal  transformation  of  the  coordinate 
axes;  it  simply  allows  us  to  use  the  experimental  values  as  a  metric  by  which  the  accuracy  of  our 
computed  elastic  constants  can  be  measured.  (We  have  computed  values  for  C14  using  the  proper 
rhombohedral  setting,  however.)  The  theoretical  and  experimental  elastic  constants  are 
presented  in  table  3.  The  agreement  with  experimental  values  is  excellent;  however,  it  should  be 
noted  that  our  theoretical  values  are  for  idealized  structures  with  perfect  B4C  stoichiometry, 
whereas  the  experimental  values  were  obtained  from  single  crystalline  measurements  of  a  sample 
with  stoichiometry  B5.6C  (6).  This  will  invariably  affect  the  experimental  results  as  compared  to 
a  perfect  B4C  stoichiometry,  though  it  is  not  clear  to  what  extent. 
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The  pressure  evolution  of  the  elastic  constants  over  the  0-80  GPa  range  for  each  polytype  is 
shown  in  figure  3.  There  we  show,  for  the  first  time,  that  for  the  BC  ceramic  each  Cij  increases 
with  pressure  however  the  C44  and  C66  moduli  show  a  marked  decrease.  This  pressure  softening 
of  elastic  moduli  (C44  in  particular)  has  been  observed  experimentally  (as  well  as  theoretically)  in 
alpha  quartz  which  is  known  to  undergo  pressure  induced  amorphization  similar  to  the 
phenomenon  observed  in  BC. 

Table  3.  Theoretical  and  experimental  elastic  constants  (GPa). 

(The  non-zero  Ch  modulus  is  indeterminate  in  this  orientation.) 


Cij 

ccc 

Polar 

Equatorial 

Exp. 

Cll 

487 

554 

559 

543 

C12 

117 

121 

117 

131 

C13 

66 

65 

70 

64 

C33 

525 

526 

521 

535 

C44 

133 

155 

153 

165 

C66 

183 

216 

218 

206 

6 


Equatorial 


Pressure  (GPa) 


Polar 


Figure  3.  Pressure  evolution  of  the  elastic  moduli  for  each  polytype. 
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The  lowest  eigenvalue  of  the  elastic  constant  tensor  for  each  polytype  is  plotted  as  a  function  of 
pressure  in  figure  4.  For  each  polytype  the  eigenvalue  decreases  monotonically  towards  zero  and 
reaches  zero  between  44-46  GPa  for  the  C-C-C  polytype  and  between  66-68  and  68-70  GPa  for 
the  equatorial  and  polar  polytypes,  respectively,  indicating  that  the  C-C-C  polytype  reaches  an 
elastic  instability  at  much  lower  pressures  than  the  other  polytypes  and  is  the  initial  point  of 
failure  in  the  BC  ceramic  under  hydrostatic  loading.  This  is  in  support  of  the  work  of  Fanchini 
et  ah,  who  studied  the  stability  of  BC  polytypes  relative  to  the  segregated  boron  and  carbon 
phases  (77).  In  their  work  they  demonstrated  that  the  C-C-C  polytype  is  relatively  unstable 
energetically  as  compared  to  the  polar  and  equatorial  configurations  as  a  function  of  hydrostatic 
load. 


Lowest  Eigenvalue 


—  ccc 

—  Polar 
Equatorial 


Figure  4.  Lowest  eigenvalue  of  elastic  constant  tensor  for  each  polytype. 


4.  Conclusions 


Our  calculations  thus  far  have  outlined  the  mechanical  response  of  BC  under  hydrostatic 
compression  as  a  function  of  polytype;  however,  it  has  been  suggested  that  shear  plays  a  critical 
role  in  weakening  the  BC  ceramic.  Our  preliminary  results  lend  support  to  this  hypothesis. 
Specifically,  examination  of  the  eigenvectors  associated  with  the  soft  mode  eigenvalues  plotted 
in  figure  4  show  a  significant  contribution  from  shear  and  indicates  that  loading  along  these 
shearing  modes  may  hasten  amorphization.  Future  work  will  attempt  to  construct  an 
amorphization  contour  in  a  selected  two-dimensional  (2-D)  pressure — shear  space.  It  now 
appears  feasible  to  use  theoretical  constitutive  representations  from  nonlinear  elasticity  to 
minimize  the  QM  calculations  so  that  full  exploration  of  the  six-dimensional  (6-D)  stress  space  is 
unnecessary. 
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The  stability  analysis  conducted  herein  allows  for  identification  of  the  stress  at  which 
amorphization  may  initiate;  however,  one  of  our  ultimate  goals  is  to  relate  changes  in  the  local 
atomic  structure  to  the  root  cause  of  the  instability.  For  example,  Binggeli  et  al.  (7)  studied  the 
pressure-induced  amorphization  phenomenon  in  a-quartz  and  presented  evidence  that  a  shear 
instability  was  the  driving  force  toward  collapse  of  the  material.  Specifically,  when  subjected  to 
shear,  silicon  atoms  were  displaced  towards  oxygen  atoms,  leading  to  a  change  in  silicon 
coordination.  Previous  theoretical  work  on  BC  showed  that  under  uniaxial  loading  the 
icosahedra  remain  relatively  unaltered;  however,  there  was  a  bending  of  the  3-atom  chain  that 
eventually  led  to  the  formation  of  new  covalent  bonds  between  the  chain  atoms  and  boron  atoms 
in  the  icosahedra  (18).  This  type  of  structural  rearrangement  can  lead  to  an  irreversible  alteration 
of  the  lattice,  which  may  trigger  the  collapse  of  BC  leading  to  amorphization.  This 
rearrangement  mechanism  needs  further  study  and  molecular  dynamics  simulation,  where  the 
change  in  structural  morphology  can  be  studied  in  real  time,  is  an  ideal  method  for  this  analysis. 
However,  this  is  entirely  contingent  upon  the  availability  of  an  accurate  interatomic  potential 
applicable  to  BC. 

This  work  has  shown  that  the  C-C-C  polytype  is  in  fact  the  weakest  polytype  within  the  BC 
structure;  however,  these  results  are  based  on  idealized  structures  with  perfect  B4C 
stoichiometry.  In  the  actual  material,  there  is  a  high  degree  of  substitutional  disorder  as  well  as 
defect  sites  consisting  of  vacancies  in  the  atomic  chain  linking  the  icosahedra.  A  computational 
survey  of  the  possible  defect  polytypes  within  BC  has  been  conducted  by  Saal  et  al.  (19); 
however,  their  worked  focused  on  structure,  enthalpy  of  fonnation,  and  phonon  spectra. 
Therefore,  calculations  of  the  mechanical  properties  of  these  additional  structures  are  required 
for  a  complete  stability  analysis  of  the  BC  ceramic. 
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Advanced  Ceramics  and  Composites  (Taylor,  201 1)  and  at  the  American  Society  of  Mechanical 
Engineers  (ASME)  Applied  Mechanics  and  Materials  Conference  (Wright,  2011). 
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